In [1]:

    

%matplotlib inline

import ephem

import numpy as np
import matplotlib.pyplot as plt

.95140571


    

    
Out[1]:





0.95140571

In [2]:

    

tles = """0 ES'HAIL 2
1 43700U 18090A   18319.95140571 -.00000497  00000-0  00000+0 0  9996
2 43700  25.0155 199.6281 7401796 178.1843  52.4356  2.15442580    02
"""

sat = ephem.readtle(*tles.split('\n')[:3])


    

In [3]:

    

mu = 3.986004418e14

def epochmjd(sat):
    return sat._epoch - 14980 

def sma(sat):
    n = sat._n / (3600 * 24) # in revolutions / second
    a = (mu/(4*np.pi**2*n**2))**(1/3) # in metres
    return a * 1e-3 # in km

def ecc(sat):
    return sat._e

def inc(sat):
    return np.rad2deg(sat._inc)

def raan(sat):
    return np.rad2deg(sat._raan)

def aop(sat):
    return np.rad2deg(sat._ap)

def ma2ta(ma_deg, eccentricity, kepler_iterations = 100):
    M = float(ma_deg)
    e = eccentricity
    # Solve Kepler equation for eccentric anomaly by iteration
    E = M
    for _ in range(kepler_iterations):
        E = M + e * np.sin(E)
    nu = 2*np.angle(np.sqrt(1-e)*np.cos(E/2) + 1j*np.sqrt(1+e)*np.sin(E/2))
    return np.rad2deg(nu)

def ta(sat):
    return ma2ta(sat._M, sat._e)


    

In [4]:

    

def gmat_keplerian_spacecraft(sat, spacecraft_name, mass = 3000, drag_area = 15, srp_area = 15):
    name = spacecraft_name
    return f"""Create Spacecraft {name};
{name}.DateFormat = UTCModJulian;
{name}.Epoch = '{epochmjd(sat)}';
{name}.CoordinateSystem = EarthMJ2000Eq;
{name}.DisplayStateType = Keplerian;
{name}.SMA = {sma(sat)};
{name}.ECC = {ecc(sat)};
{name}.INC = {inc(sat)};
{name}.RAAN = {raan(sat)};
{name}.AOP = {aop(sat)};
{name}.TA = {ta(sat)};
{name}.DryMass = {mass};
{name}.DragArea = {drag_area};
{name}.SRPArea = {drag_area};
"""


    

In [5]:

    

print(gmat_keplerian_spacecraft(sat, 'Eshail2'))


    

    




Create Spacecraft Eshail2;
Eshail2.DateFormat = UTCModJulian;
Eshail2.Epoch = '28438.451405710002';
Eshail2.CoordinateSystem = EarthMJ2000Eq;
Eshail2.DisplayStateType = Keplerian;
Eshail2.SMA = 25322.941234489062;
Eshail2.ECC = 0.7401795983314514;
Eshail2.INC = 25.015499114990234;
Eshail2.RAAN = 199.6280975341797;
Eshail2.AOP = 178.18429565429688;
Eshail2.TA = 140.81615297409934;
Eshail2.DryMass = 3000;
Eshail2.DragArea = 15;
Eshail2.SRPArea = 15;

In [6]:

    

sma(sat) * (1 + ecc(sat))


    

    
Out[6]:





44066.46570600413

In [7]:

    

earth_radius = 6371
sma(sat) * (1 + ecc(sat)) - earth_radius


    

    
Out[7]:





37695.46570600413

In [8]:

    

sma(sat) * (1 - ecc(sat)) - earth_radius


    

    
Out[8]:





208.41676297399954